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1  Introduction 

Aerodynamic  shape  design  is  a  challenging  application  for  modem  computational  fluid  dynamic.  This  is,  in¬ 
deed,  a  very  productive  field,  both  in  term  of  scientific  publications  and  developed  applications.  Furthermore 
it  is  possible  to  find  excellent  analysis/design  software  released  in  the  public  domain  (See  for  example  [1,  2]). 
However,  industry,  in  particular  automotive  and  aerospace,  but  also  emerging  fields  like  wind  and  sea  extracted 
energy,  continuously  needs  increased  product  performance  and  competitiveness.  This  requires  increasingly  re¬ 
fined,  advanced  and  complex  computational  models  capable  of  describing  flow  behavior  with  higher  and  higher 
fidelity. 

Thus  the  design  tools  available  to  scientist  and  engineers  arc  challenged  by  the  growing  complexity  of  sim¬ 
ulation  models.  The  increasing  price  to  performance  ratio  of  nowadays  computers  is  a  strong  allied  of  designers 
in  this  though  race,  but  often  it  is  not  enough. 

Objectives  of  this  lecture  is  to  explore  and  analyze  some  of  the  technologies  that,  together  with  the  increas¬ 
ingly  available  computational  power,  may  help  to  exploit  the  tremendous  potential  of  high  fidelity  flow  field 
models  to  aerodynamic  shape  design  challenges. 

For  the  sake  of  simplicity  we  will  concentrate  on  aerofoil  design.  Two  aspects  will  be  considered:  the  com¬ 
puter  algebra  assisted  development  of  code  related  to  gradient  computation  through  adjoint  techniques  and  the 
use  of  fitness  approximation  techniques  to  improve  the  performance  of  search  procedures.  The  concern  of  this 
first  lecture  will  be  more  methodological,  while  some  applications  will  be  illustrated  in  the  next  one. 

The  outline  of  the  first  first  part  of  this  lecture  is  as  follows:  some  short  remarks  on  adjoint  problems  will 
introduce  the  framework  in  which  the  Computer  Algebra  techniques  will  be  mainly  used  here.  Then  a  brief 
outline  of  symbolic  computation  systems  will  be  given,  with  particular  reference  to  the  features  that  arc  here  used 
for  adjoint  equation  derivation.  The  implemented  algorithm  will  be  described  and  illustrated  with  an  example 
on  a  simple  problem.  The  adjoint  and  the  related  boundary  conditions  for  3D  Navier-Stokes  arc  then  reported, 
while  an  example  of  derivation  is  given  in  appendix  A. 

The  second  paid  of  this  lecture  will  focus  on  the  use  of  approximated  fitness  evaluators.  The  problem  will 
be  analyzed  using  the  framework  of  multi-objective  optimization,  and  some  possible  approaches  to  the  problem 
of  mixing  exact  and  approximated  evaluations  will  be  discussed. 


2  Adjoint  formulation  for  gradient  evaluation 

The  adjoint  formulation  appeal's  in  optimal  control  theory  as  a  tool  to  compute  the  constrained  gradient  to  a 
given  functional.  In  fluid  mechanics  the  applications  span  from  shape  design  to  flow  control  [3,  4,  5,  6].  The 
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adjoint  method  is  based  on  the  solution  of  an  additional  set  of  equations  for  the  Lagrange  multipliers.  The  La¬ 
grange  multipliers  local  value  is  a  measure  of  the  functional  sensitivity  to  the  flow  variables  local  variation.  This 
provides,  for  example,  a  quantitative  criteria  for  grid  refinement  [7], 

In  mathematical  terms  the  adjoint  approach  to  gradient  computation  can  be  concisely  derived  from  the  La¬ 
grange  identity: 

(. Av,w )  =  ( v,A*w )  (1) 

where  (  ,  )  is  the  scalar-  product  in  the  appropriate  Hilbert  space.  This  identity  defines  the  adjoint  operator 

A*  when  a  linear-  or  non-linear-  operator  A  and  when  v  £  D,  w  G  D*  are  given.  Now,  if 

J  =  (p,v)  (2) 


is  a  given  functional  and 


Av  =  /,  A*w  =  p 


(3) 


are,  respectively,  the  state  and  the  adjoint  equations,  then  it  is  possible  to  write 

J  =  (p,  v)  =  (/,  w)  (4) 

as  a  direct  consequence  of  identity  1 .  This  means  that  J  can  be  expressed  either  in  terms  of  the  state  vector  v  or 
in  terms  of  the  adjoint  one  w.  The  small  perturbations  theory  can  then  be  used  to  obtain  the  functional  variation 
for  non-linear-  state  equations.  In  general,  if  A(v)  is  an  operator  that  depends  non-linearly  on  v,  it  is  possible  to 
write,  under  the  hypothesis  of  small  perturbations,  that 

/  =  /o  +  £/i  +  •  •  • 


v  =  v0  +  evi  +  . . . 

d  A 

A(v0  +  evi  +  . . . )  =  2f(u0)  +  e  — 

av 


Vl  +  . . . 


Therefore  the  non-linear-  state  equation  A(v)v  =  f  can  be  decomposed  as  follows 

-4(ToH  =  /o 


vi  Wo  =  fl 


A(v,)Vl  +  ^ 

From  the  last  equation  descends  immediately  the  first  order  linearized  operator: 

dA\ 


-4i(^o)«  =  A(vq)  •  + 


dv 


vo 


In  the  same  way,  it  is  possible  to  write  the  functional  J  as 

J  =  Jo  +  £</i 

where 

Jo  =  (vo,p)  =  (wo,fo) 
§J  =  J1  =  (vi,p)  = 


(5) 

(6) 

(7) 

(8) 

(9) 


with  AI(vq)wi  =  p 

In  other  words  we  obtain  explicitly  the  functional  variation  with  respect  to  the  right  hand  side  (RHS)  of  the 
governing  equation,  i.e.,  the  gradient  of  J  with  respect  to  /.  Hence,  from  a  mathematical  point  of  view,  the 
derivation  of  the  adjoint  state  equations  for  a  given  functional  J  is  a  straightforward  process.  A  disadvantage  of 
this  method  is  that  complicated  functionals  often  require  a  lengthy  and  error-prone  manual  derivation  and  hand 
coding  process.  In  this  respect  Computer  Algebra  (CA)  techniques  can  be  profitably  used  to  ease  and  speed-up 
this  process.  In  the  following  sections  such  technique  will  be  described  and  then  used  to  automatically  obtain 
the  adjoint  equations  and  pertinent  boundary  conditions  of  the  Navier-Stokes  equations.  Gradient  expression  is 
also  automatically  computed  for  a  generic  objective  function.  A  standard  computer  algebra  language  (MAPLE 
[8])  is  used  to  build  the  automated  derivation  procedure. 
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3  Automatic  derivation  of  the  adjoint  system  and  gradient 

3.1  General  remarks  on  symbolic  computation 

Computer  algebra  systems  for  symbolic  computation  have  a  number  of  peculiarities  that  make  them  radically 
different  from  numeric  computing  tools.  Indeed,  they  operate  on  symbols  (the  rich  set  of  mathematical  objects), 
rather  than  numbers  as  classic  numerical  programming  languages.  The  mathematical  objects  arc  transformed 
following  exact  equivalence  rules,  and  there  are  not  precision  limits  in  the  representation  of  numbers.  The  al¬ 
gebraic  manipulation  algorithms  can  be  very  complex  and  sophisticated  so  that  a  computer  algebra  language 
can  be  considered  a  way  of  storing  and  efficiently  using  the  knowledge  and  expertise  of  thousands  of  brilliant 
mathematicians.  This  allows  an  easy  and  straightforward  use  of  computer  algebra  systems  as  “scratchpad”  tools 
in  the  interactive  solutions  of  complex  mathematical  problems. 

On  the  other  hand,  there  arc  some  drawbacks  that  have  to  be  taken  into  account  when  developing  complex 
applications.  In  particular,  there  arc  cases  in  which  the  output  is  very  complex,  and  considerable  effort  has  to 
be  devoted  to  keeping  it  easily  understandable.  Furthermore,  it  has  to  be  also  considered  that  it  is  not  easy  to 
learn  how  to  write  efficient  code  in  a  computer  algebra  language,  because  the  programming  model  is  fundamen¬ 
tally  different  with  respect  to  “numerical  programming”,  with  recursion,  list  and  set  operations,  and  expression 
identification  and  and  matching  on  a  complex  and  rich  set  of  objects. 

The  approach  followed  in  the  present  work  was  the  research  of  a  compromise  between  output  conciseness 
and  ease  of  implementation,  and  the  main  steps  of  the  developed  adjoint  equation  derivation  procedure  can  be 
summarized  as  follows: 

•  The  Navier-Stokes  equations  have  been  specified  term  by  term,  but  the  volume  and  surface  integrals  of 
the  Lagrangian  have  been  left  unexpanded. 

•  The  classical  manipulations  for  adjoint  derivation  were  obtained  using  Green  formulas,  while  the  variation 
with  respect  to  geometry  were  obtained  through  Hadamard  formulas. 

This  approach  required  the  capability  of  discriminating  between  the  various  possible  structures  of  a  given  La¬ 
grangian,  in  order  to  apply  the  proper  transformations.  This  was  obtained  applying  recursively  the  pattern  match¬ 
ing  functions  of  the  CA  system. 

3.2  Adjoint  computation 

We  provide  here  expressions  and  identities  that  will  be  useful  for  the  automatic  computation  of  the  integral  vari¬ 
ations  by  the  MAPLE  V  computer  algebra  system  [8]. 

The  first  step  in  the  automated  derivation  process  is  the  expression  of  the  flow  field  equation  in  indexed  form, 
using,  in  the  present  implementation,  the  Grtensor  software  [9]  for  tensor  calculus.  The  expanded  equation  and 
boundary  conditions,  referenced  as  R,  are  then  combined  with  the  objective  function  O  to  form  the  Lagrangian 
equation 


L  =  1)0(  X ,  •  •  •  ,  n,  Uxi  'Uyi  'U’Zi  ^ xxi  ^ xyi  •  •  •  i  •  •  •  )  T 
£  (h) 

f 

AR(x, ...  ,  rr,  Ux^  Uzi  ^ xxi  ^ xyi  •  •  •  ,  h: . . . )  dfl 


n  (h) 


(10) 


The  adjoint  equations  and  boundary  conditions  arc  obtained  computing  the  variation  of  L  with  respect  to  the  flow 
field  unknowns  u,  while  the  gradient  has  been  computed  through  the  variation  of  L  with  respect  to  the  boundary 
defining  functions  h.  The  formulas  used  for  the  variation  of  volume  and  surface  integrals  arc  reported  in  the 
following  sections.  One  of  the  tricky  aspects  of  this  approach  is  that  the  variation  of  L  has  to  be  computed  with 
respect  to  quantities  that  arc,  in  turn,  functions  of  the  coordinates  x,  y  and  z.  In  fact,  standard  MAPLE  differ¬ 
entiation  function  is  unable  to  perform  variational  calculus.  Yet,  there  arc  ways  to  get  around  such  a  difficulty. 
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for  example  using  the  pdif  f  [10]  function,  available  in  MAPLE  V  share  library.  Such  function  enables  the 
computation  of  expressions  such  as 

df{x,u{x )) 
du{x) 

that  arc  not  computable  with  the  standard  MAPLE  differentiation  function.  The  Green’s  theorem  is  then  repeat¬ 
edly  applied  to  move  the  terms  containing  the  derivative  of  a  variation  from  the  volume  terms  to  the  surface 
integral  terms.  MAPLE  pattern  matching  functions  have  been  extensively  used  both  for  the  variation  procedure 
and  the  Green’s  theorem  procedure  implementation. 

3.3  Symbolic  computation  algorithm 

Given  the  Lagrangian  expression,  in  the  following  we  give  some  details  about  the  algorithm  used  to  obtain  the 
adjoint  equations,  the  boundary  conditions,  and  the  gradient  expression.  The  derivation  of  the  formulas  adopted 
to  express  the  variation  of  the  functionals  is  reported  in  appendix  B. 

The  first  MAPLE  function  implemented  is  called  variation.  This  function  takes  as  input  a  generic  func¬ 
tional  L,  a  function  u  to  be  variated  and  eventually  the  vector  function  h  that  defines  the  geometric  deformation. 
In  general,  the  expansion  obtainable  for  a  generic  variation  is  of  the  kind 

8Lu(x,  .  .  .  ,  U,  Ux,  Uyj  Uzj  'U’xxi  ^xyi  •  •  •  )/&)•••)  — 


dL  .  dL  . 

—  drt-l-  - — dux 
du  dux 


dL  A 

d^Suy 


dL 

duz 


.  dL  .  dL  . 

OUz  +  — - OUxx  +  — - OU: 


du 


du  ~“'xy 

L/  UXy 


(11) 


When  the  variation  is  computed  with  respect  to  the  geometry,  a  pattern  matching  algorithm  checks  each  term  of 
the  functional  L  in  order  to  find  volume  or  surface  integrals.  In  the  case  of  a  surface  integral,  then  the  following 
expression  is  used  to  compute  the  variation: 

SF  =  JJ  V/-^<5pdrr  +  II  fU~8gda  +  H^dgda  (12) 

s  s  s 


Appendix  sub-section  B.l  reports  the  derivation  of  this  expression  and  the  meaning  of  the  terms  involved.  In 
the  case  of  volume  integrals,  the  formula 


"-JIJ  ffW/'S-** 

T  FT 


(13) 


is  adopted.  Expression  derivation  and  terms  definition  arc  in  B.3.  These  formulas  arc  exactly  equivalent  to  the 
Hadamard  formulas.  It  is  worth  to  remark  here  that  volume  and  surface  integrals,  need  special  treatment  only 
when  the  integration  manifold  is  given  in  implicit  form.  In  the  other  cases  they  can  be  reduced  to  equivalent 
double  and  triple  definite  integrals  and  treated  with  formula  1 1 . 

The  second  MAPLE  function  implemented  is  called  intparts3d.  The  input  to  this  function  is  the  vari¬ 
ated  functional  5L,  and  in  output  an  equivalent  form  of  the  functional  5L  is  generated,  were  the  terms  contain¬ 
ing  derivatives  of  the  variation  arc  moved  to  surface  integrals.  It  works  applying  recursively  the  integration  by 
parts  and  the  Green’s  theorem  to  the  volume  integral  terms  containing  derivatives  of  the  variation,  e.g.  8ux  or 
8uxy.  The  resulting  functional  variation  presents  terms  containing  derivatives  of  the  variation  only  inside  sur¬ 
face  integrals.  In  a  subsequent  step,  these  terms  can  be  eliminated  through  the  imposition  of  suitable  boundary 
conditions.  This  last  step,  currently,  is  not  automated  and  require  some  care  and  attention  on  the  pari  of  the 
user.  To  help  this  process,  the  functions  extract_ad  j_eq  and  extract_ad  j  _bceq  have  been  introduced. 
They  extract  the  adjoint  equations  and  boundary  conditions,  respectively,  again  by  pattern  recognition.  Finally, 
two  interface  procedures,  compute_ad  j_eq  and  compute_ad  j_eq_bc,  hevebeen  introduced  to  manage  the 
whole  process  of  adjoint  equation  and  boundary  conditions  derivation.  The  listing  in  MAPLE  code  of  the  first 
one  is  reported  below. 
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compute_ad j_eq  := 

proc (lagrangian,  v,  dv,  n_intp,  h,  verbose_level ,  aliasing :: boolean) 

option  remember; 

local  lg_v,  i,  cl,  sp,  adj_eq; 

lg_v  :=  variation ( lagrangian,  [v=v+dv] ,  h,  [si,  s2],  [xl,  x2,  x3]); 

lg_v  :=  expand ( lg_v) ; 

if  aliasing  then  autoAlias ( lg_v)  fi; 

if  verbose_level  >=  2  then  print ( "LAGRANGIAN" ) ;  print (lg_v)  fi; 
cl  : =  1 g_v ; 

for  i  from  1  by  1  to  n_intp  do 
cl  :=  intparts_3d (cl , dv) ; 
cl  :=  expand ( convert (cl ,  diff)  )  ; 
if  aliasing  then  autoAlias (cl )  fi; 
if  verbose_level  >=  1  then 

print ("INTEGRATION  BY  PARTS  N.",  i);  print (cl) 

fi; 

od; 

adj_eq  :=  extract_adj_eq (cl, dv)  ; 

end : 


The  equations  obtained  by  the  automated  procedure  arc  of  course  less  concise  compared  to  those  obtained 
manually.  Indeed,  the  governing  equations  arc  fully  expanded  as  a  first  step  in  the  variation  computation,  to 
ease  the  pattern  matching  process.  Furthermore,  the  comparison,  useful  for  debug  purposes,  with  the  manually 
obtained  equations  is  complicated  by  some  asymmetries  that  may  be  found  in  the  automatically  derived  terms. 
This  is  due  to  the  following  identity  (see  sub-appendix  B.4  for  details): 


dfl 


vxri2  dTi 


vyrii  dTi 


n 


E 


E 


(14) 


This  comes  from  the  application  of  the  divergence  theorem  to  terms  containing  crossed  derivatives.  Indeed,  the 
theorem  can  be  applied  here  in  two  different  ways,  and  the  automatic  procedure  does  not  know  which  form  has 
to  be  chosen  to  keep  symmetry.  Flowever,  the  same  identity  14  can  be  applied  to  suitably  transform  the  obtained 
expressions. 


3.4  A  simple  application  example  to  potential  flows 

An  example  of  adjoint  equation  derivation  is  here  given  to  show  how  the  developed  MAPLE  code  may  be  used. 
The  reader  may  refer  to  appendix  A  for  the  application  to  Navier-Stokes  equations. 

A  symmetric  airfoil  is  immersed  in  a  potential  flow  at  zero  angle  of  attack  (see  figure  1),  and  a  symmetric 
transpiration  velocity  W  is  assigned  in  the  direction  normal  to  the  airfoil  surface.  The  governing  equations  arc 
therefore: 


V20  =  0  in  the  field,  V</>-n  =  W  on  the  airfoil 


(15) 


The  problem  is  to  find  a  transpiration  velocity  distribution  W*  such  that  a  given  target  potential  '!>  is  matched. 
Flence  the  objective  function  that  will  be  minimized  is: 


1 

2 


/(</>  —  <b)2der 


(16) 


on  the  airfoil  surface.  The  current  version  of  adjoint  derivation  code  was  developed  taking  into  account  three- 
dimensional  fields;  therefore  in  the  following  example,  the  two-dimensional  field  will  be  treated  as  3D  field  with 
zero  gradient  of  the  flow  variables  in  x3  direction.  In  the  same  spirit  the  airfoil  is  a  cylindric  surface  identified 
by  two  curvilinear-  coordinates  si,  s2.  Thus  we  will  have  surface  integrals  instead  of  contour  integrals  on  the 
boundary  and  volume  integrals  in  the  field. 
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W  =  V0-n 


Here  the  MAPLE  interactive  session  for  the  adjoint  computation  follows.  The  alia  s_depe  nde  n  c  i  e  s 
function  set  a  functional  dependence  for  a  given  list  of  variables;  dW  is  the  increment  of  the  blowing  ve¬ 
locity,  y  is  the  airfoil  ordinate,  xl,  x2,  x3  are  the  field  coordinates,  while  si,  s2  are  the  surface  coordinates. 
Note  that  to  have  the  procedure  working  properly,  we  had  to  specify  explicitly  the  dependence  of  <j>  and  of 
its  increment  dcj)  on  field  and  surface  coordinates.  The  same  dependence  is  specified  for  the  vector  normal 
to  the  airfoil  surface  n  =  [nl,  n2,  ra3]  and  for  the  deformation  field  h  =  [hi,  h2,  h3\  (not  explicitly  used 
here). 

>  al±as_dependencies ( [y,  W, dW]  ,  [ si ] )  : 

>  varlist  :=  [xl, x2, x3, si, s2]  : 


> 

> 

> 


> 

> 

> 

> 

> 

> 

> 


alias_dependencies  (  [phi ,  dphi,  nl,n2,n3,hl,h2,h3]  ,  varlist )  : 
alias_dependencies ( [Phi]  ,  [si,  s 2 ]  )  : 

Obj :=Int (Int (1/2* (phi-Phi) ~2,sl),s2); 


Obj  :  = 


(4>  —  <f>)2  dsl  ds2 


Potential_eq : =dif f (phi,  '$'(xl,2))+diff(phi,  ' $ ' (x2 , 2 ) )  : autoAlias (% ) ; 

4^x1 ,  xl  T  4>x2,  x2 


n : = [nl , n2 ] ; 


n  :=  [nl ,  n2\ 

Neumann_bc : =dotprod (grad (phi, [xl , x2 ] ) , n, 'orthogonal' ) -W: 
autoAlias (%) ; 


<t>x!  nl  +  (j)X2  n2  -  W 


alias_dependencies ( [eta, beta]  ,  varlist)  : 

Lagrangian : =0b j+Int (Int (Int (eta*Potential_eq, xl) , x2) , x3) + 
Int (Int (beta*Neumann_bc, si ) , s2 )  ; 


Lagrangian  := 


(4>  —  <b)2  dsl  ds2  + 


V  {.4>xi ,  xi  +  4>x2  ,  x2 )  dxl  dx2  dx3 


+  /3  ( 4>xi  nl  +  4>x2  n2  —  W)  dsl  ds2 


>  ad j_dphi : =compute_ad j_eq (Lagrangian, phi, dphi, 2, [hl,h2,h3],0,true); 

adj -dphi  . —  rjxi,xi  T  iJX2,x2 

>  ad j_dphi_bc : =compute_adj_bc (Lagrangian,  phi,  dphi,  2,  [nl,  n2,  n3]  , 

>  [hi,  h2,  h3],  0,  true) ; 


adj _dphi_bc  :  = 

,  ^  y  dphir1  nl 

</>-$+'  xl 


Vxi  nl 


y  dphix2  n2 


-  yx2  n2  + 


/ 3  nl  dphixl  (3  n2  dphix2 


dphi  dphi  dphi 

ad j_dphi_bc_airf : =subs (beta=-eta, ad j_dphi_bc) ; 


dphi 
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adj -dphi-bc-airf  :=  p  —  $  —  rjxi  nl  —  rjX2  n2 

>  c6 :=variat±on (Lagrangian,  [W=W+dW]  ,  [sl,y,s2],  [sl,s2],  [xl, x2, x3] )  ; 


P  dW  dsl  ds2 


In  the  above  equation  the  gradient  with  respect  to  the  blowing  velocity  W  can  be  recognized  as  being  given 
by  the  variation  of  the  Lagrangian  with  respect  to  W. 


4  The  compressible  Navier  Stokes  Adjoint  Equations 

In  this  section  the  mathematical  adjoint  to  the  Navier-Stokes  equations  are  derived.  The  notation  is  slightly 
different  from  that  of  the  previous  chapter  for  the  sake  of  conciseness  and  clarity.  In  particular,  the  variations 
will  be  indicated  with  a  circumflex  accent. 


4.1  State  equations 

In  this  section  we  make  the  notation  clear-  by  giving  the  steady  compressible  Navier-Stokes  equations  under 
divergence  form  and  indicial  notation. 


Conservation  of  mass 


Conservation  of  momentum 


dpvj 

dxi 


=  0 


d 

dxj 


|  P'OtVj  +pdij  -  p 


dvi  dvj 
dxj  dxi 


2  dvk  A  ' 

3  dxk 


=  0 


Conservation  of  energy 


_d_ 

dxi 


n  tdT 
pviH  -  k—  -  p 


dvj 

dXj 


9vj 

dxi 


2  dvk  A  ' 

3  dxk 


=  0 


(17) 


(18) 


(19) 


H  is  the  total  enthalpy  per  unit  mass,  k  the  heat  conductivity,  and  the  other  symbols  have  their  usual  meaning. 
Let  Q  be  the  flow  domain  and  E  its  surface.  The  no-slip  boundary  conditions  are  imposed  on  E,  elsewhere 
we  consider  suitable  far-  field  conditions.  The  boundary  condition  on  temperature  is  to  be  detailed  later. 


4.2  Lagrangian 

I  is  a  generic  functional  depending  on  the  flow  field  and  it  is  defined  either  over  Q  or  over  E.  We  want  to  find 
a  surface  E  such  that  the  first  order  variation  of  /  with  respect  to  a  perturbation  of  E  is  zero. 

To  this  end,  we  introduce  an  auxiliary  functional  L  defined  as  follows 


L  =  I  + 


+ 

+ 

+ 

+ 


f  dpv, 

LV  dxi 
d 


1  dn 


in  dxn 

f  c— 

In  S  dxi 


A*  —  yrojVj  +  p5i  j  -  p 

IT  l9T 
pviH  -  k—  -  p 


dvj  dvj_  _  2  dvk 

dxj  dxi  3  dxk 

dvj  dvj_  _  2dvk_  . 

dxj  dxi  3  dxk 


Vj  >  dfl 


dE 

t0(T)  dS 


(20) 
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0(T)  =  T  —  Tw  with  prescribed  boundary  temperature,  whereas  for  adiabatic  flows  ©(T)  =  dT/dxiUi. 
The  auxiliary  functions  77,  A*  and  £  belong  to  H2{ fl).  £  £  L2(S)  and  r  £  L2(S)  or  r  £  Hl(T)  according  to 
the  boundary  conditions  selected  for  temperature.  We  refer  to  all  these  functions  as  Lagrange  multipliers. 

The  lagrangian  function  L  transforms  the  constrained  problem  defined  for  I  into  an  unconstrained  problem: 
it  is  now  required  to  find  the  flow  field,  E  and  the  Lagrange  multipliers  in  order  that  the  respective  first  order 
variations  of  L  arc  zero.  The  variation  of  L  with  respect  to  the  Lagrange  multipliers  yields  the  flow  equations 
and  boundary  conditions.  We  concentrate  on  the  variation  of  L  with  respect  to  the  conservative  flow  variables 
p,  pvi,  pE,  with  E  the  total  energy  per  unit  mass. 

Take  the  variation  of  the  first  integral  in  eq.  20,  using  the  Gauss  identity  we  have 


p  pvi  rii  dT  — 


In 


~  dp 
pvi  — —  di  l 
OXi 


(21) 


where  rii  is  the  outward  unit  normal  to  E,  and  •  denotes  the  function  variation. 
Consider  now  the  second  integral  in  eq.  20  and  apply  the  Gauss  identity  twice  to  get 


+ 


in 


(, pviVj  +  pSi  j)  dQ  - 


dX3 

dxi 


2d\ k 
3  dxk  12 


ViUj  dT, 


dxi 


2  d\k  s 

3  dxk 


Vi  dfl 


(22) 


where  the  components  of  the  force  per  unit  surface  experienced  at  the  boundary  arc 


ft 


|  pvp)3  +p5ij  -  p 


dvi  dvj 
dxj  dxi 


2  dvk  x 

3  dxk 


n, 


To  be  remarked  is  that  the  linear-  second  order  adjoint  operator  keeps  the  same  form  of  the  viscous  stresses. 
In  eq.  22  the  flow  variables  do  not  appear-  under  the  form  of  conservative  variables.  However,  in  the  first 
integral,  for  a  reason  that  will  be  clear-  later,  we  leave  the  flow  variation  under  the  form  /,,  while  in  the  others 
we  write  the  variations  in  terms  of  conservative  variables.  For  example  we  have 


pViVj  =  pViVj  +  pVjVi  -  ViVjP 


1  Vi 

Vi  =  -pvi - p 

p  p 


p  = 


7-1 


(2pE-2  pvm  +  V2p) 


T  = 


7-1 
2  pR 


2  pE 


2  pviVi  + 


P 


pviH 


pviH  +  7 VipE  -  (7  -  1  )vipvjVj  +  Vi 


P 


with  7  the  specific  heats  ratio,  R  the  perfect  gas  constant  and  V2  =  7777.  Substituting  the  above  relations  in 
eq.  22  and  posing 


hi 


d  (  d\j 
dxj  \  dxj 


d^j 

dxi 


2d\k  1\ 

3  dxk  tJ\  J 
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one  finds 


A  ifidT- 

7-1 


Pi 

dX , 

„  -VjVj  — — 
2  J  3  dxi 


dX i  d\j  .  dXj  hi 

v*87j+v*d£-{',-1)v‘ed  +  7. 


pvi  dXl 


dXj 


1 


-  ViVj— - -Vihi  )  pdQ 


dxi 


a\ _ 

(7-1  )^pEdn  + 

>  UXj 


p 


dX  i  dXj 
dxj  dxi 


2  dXk 

3  dxk 


The  third  integral  in  eq.  20  is  rearranged  as  above  to  get 

3C 


f  dT 

c edT-  /  k(-dT- 
.It  oxi 


'n 


pviH-^-dQ  — 

OXj 


kT^-m  dT 

OXj 


ViTij  dT, 


f  d 

L  T~dxi 


(23) 


k  )  dn 

OXj 


+ 


P 


K 

dxj 


dvi  dvj 
dxj  dxi 


2  dvk  A 
3dxk  ij 


Vj  dfl 


(24) 


with 


e  =  pviH  -  p 


dvi  dvj 
dxj  dxi 


2  dvk 

3  dxk 


Vi  >  m 


The  variation  of  eq.  24  is  complicated  by  the  presence  of  non-linear  terms  involving  both  the  Vj  and  the 
derivatives.  As  a  first  step,  we  write  the  variation  of  eq.  24  as 


C edT —  f  k(— — riidT—  f  pv{H 

OXi  Jn 


d( 


+ 


+ 


P 


<K_ 

dxj 


pm 


d( 


dxj 


'£ 
dv, 

dxi  dx 

<K_ 

dxj 


i  dvi 
+  J 


Vi  +  —Vj 


2  dv±r 

3  dxk 


dxi 

_d_ 

dx 


dXl 


kT^mdT-  [  j 
dxi  Jq 

d(  d(  2 

^v’d^+v‘a^-3vk 


dxi 

0C*.. 

dxk  lj 


k  )  dn 

OXi 


Vjd£l 


2  d( 


v k  dj  j  Vj  dT 


(25) 


3  dxk 

Leaving  unchanged  the  terms  on  the  boundary  and  expressing  the  other  variations  in  terms  of  the  conserva¬ 
tive  variables  variations,  the  equation  above  finally  becomes 


f  dT 

(edT  -  /  k(— —  dT 

Jt.  dxi 


kT^-m  dT 

dxi 


+  J 
+  J 
~  J 

where 


9jVj  _  7-1 
p  2  pR 

7  —  1  d 

pR  dxi 


V2  - 


dxi 

2  R 


7-  1 


pR  dx. 
d 


pm 


d( 


<K 

dxj 

9j 


K 

dxi 


Vi  +  —Vj  - 


2  d( 


3  dxk 


vkSi  j  Vj  dT 


dC  'y  —  1  f) 

[(7  -  1)  VjVj  -  H5ij }  Vi —  PVJ  dn 


9 3  = 


d_ 

dxi 


K 

dxi 

d( 


IVi 


K 

dxi 


d  (  2 


1 dx 


pE  d£l 
d( 


T  Vi- —  - 


oc 


dxi 


pd9l 


dxklj 


~  P 


<K_ 

dxi 


dvi  dvj 
dxi  dxj 


2  dvk  ' 
3dxk 


(26) 


Note  the  self-adjointness  of  the  linear-  part.  The  variations  of  the  last  two  integrals  in  eq.  20  are  easily  com¬ 
puted  as 


(pH  dT 


(27) 


r0  dT 


(28) 


where  ©  =  T  or  0  =  dT/dxim ■  It  is  now  completed  the  variation  of  eq.  20.  Consider  all  the  field  integrals 
and  factor  out  the  conservative  variables  variation. 
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4.3  Adjoint  equations  and  boundary  conditions 

In  order  to  be 


L  =  0  vjp,  pvi,  pE^ 

the  following  equations  arc  to  be  satisfied 


7  ~  1 
2 


d\ i  d\ i 

VjVj- - 1"  vivj  a — 

OXn  OXi 


„.  dc  ,  vi(hi  +  9j) 

ui  o  i 

OXi  p 


( ^ ~  1 
V  2p-R 


I/2 


r\  _d_ 
p  /  5®* 


(29) 


[(7  -  1)  ViVj 


dp  dX  i 
~dxi  ~  Vj~dxj  ~  Vj 


dXj 

&-  +  <7 


hi+9i  ,  7-1  d  f  J  dC  \  _  n 

~^  +  ~^fVid^  \kd^)~° 

(30) 


t  \  5A,;  _  <9£  7-15 

dxt  1VidXi  pR  dxi 


(31) 


These  arc  the  compressible  Navier-Stokes  adjoint  equations.  The  boundary  conditions  to  be  satisfied  by  such 
equations  arc  determined  by  considering  the  boundary  terms  of  L.  We  focus  on  wing-like  geometries,  such  that 
in  the  far  field  the  flow  is  considered  unperturbed,  and  the  no-slip  boundary  condition  is  applied  on  a  finite  closed 
surface.  In  this  respect  L  depends  on  the  boundary  terms  left  in  eqs.  21,  23,  26,  27,  28,  once  the  field  integrals 
arc  elided  by  satisfying  the  adjoint  equations.  We  arc  left  with 


L  =  I 

+  / 


Ai/i<E 


dT 


d( 


s 

dx 


rjpviriidTi+  /  £ edE  —  k  /  C~^ZTni  dS  +  k  /  T— — n^dE 

is  dxi 


dXj 


2dXi 


£i  +  vpni  +  p  Tr—rij  +  -rp-rij  -  n*  + 


dxj '  "J  dxi  3  3  dxj ' 


1  s  ’dxi 

5£ 

dx, UJ 


2  d£ 

( 

Suppose  we  give  a  certain  wall  temperature.  Then  0  =  T,  and  therefore  by  taking 


<K_ 

13  dxj 


6  1  dxj  3 


j  r©dS 

Ui  dE  (32) 


dC 

£  =  0  and  r  =  —  k— — rii 
dxj 


Prescribed  wall  temperature 


(33) 


on  E,  the  boundary  integrals  involving  T  =  0  and  dTdxi  arc  eliminated.  In  contrast,  the  adiabatic  wall 
corresponds  to  ©  =  (dT / dxi)nt,  and  hence  if 


Tr~ni  =  0  and  r  =  /;:£  Adiabatic  wall 

OXi 


(34) 


on  E,  the  respective  boundary  integrals  arc  eliminated  as  in  the  previous  case.  As  vt  =  0  on  the  solid  walls,  we 
have  remarkable  simplifications  for  e: 


e  = 


|  pHn%  -  p 


dvj 

dxj 


d^ 

dxi 


2  5nfc  ' 

3  dxk  1 3 


rij 


Vi 


Also  the  third  term  in  32  is  0  and  if  we  take 


& 


pHrii  -  p 


dvi 

dxi 


2  dvk  \ 

3  dxk  ij) 


(- d 


5A  j 

dxi 


2dXk  \ 
3  dxk  ijJ 


rij  -  pnpi 


(35) 
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we  are  left  with 

L  =  I  +  J^\JidX  (36) 

where  /  is  the  variation  of  the  force  per  unit  surface  acting  on  the  fluid  at  the  boundary. 

Such  variation  is  unconstrained,  and  there  arc  two  different  ways  to  elide  the  two  terms  left  in  the  equation 
above.  The  first  is  that  I  can  be  put  under  the  form 

I  =  fjJidZ  (37) 

so  that  taking 

Aj  =  -Ji  (38) 

the  variation  of  the  lagrangian  is  finally  zero.  In  turn  this  condition  shows  how  the  only  allowable  surface  cost 
functionals  arc  those  involving  /*,  as  other  functionals  will  not  lead  to  L  =  0.  For  a  deeper  discussion  see  [11]. 
Otherwise,  if  I  is  a  field  integral  then 

i=  (  J{p,m,pE}fdn  (39) 

Jn 

where  in  this  case  J  is  the  derivative  of  the  integrand  with  respect  to  the  conservative  variables.  The  three 
terms  under  integral  become  source  terms  of  the  adjoint  equations,  and  by  taking  A  j  =  0  on  E  we  are  left  again 
with  L  =  0. 


4.4  Gradient 


As  mentioned,  the  variation  of  L  with  respect  to  the  Lagrange  multipliers  yields  back  the  governing  equations 
and  boundary  conditions,  therefore  the  gradient  is  simply  the  variation  with  respect  to  the  geometry.  If  /  is  a 
field  integral,  and  if  e  cu  is  variation  in  the  direction  of  the  normal  to  the  boundary,  we  have 

VI  =  [  ^ Ludn+  [  Flu  dT.  (40) 

Jn  olu  7s 

where  F  is  the  integrand  of  I  and  e  is  small.  When  I  is  defined  on  the  boundary  we  have 

VI  =  [  ^-LudFI+  [  |^n;u;dE  +  [  HFcudX  (41) 

J s  ou  7s  oxi  7s 

and  H  is  the  local  surface  curvature. 

For  the  other  terms  of  L  the  derivation  rules  arc  the  same,  but  some  simplifications  occur  since  the  governing 
and  adjoint  equations  with  the  respective  boundary  conditions  arc  satisfied  on  E.  The  only  terms  left  arc  those 
relative  to  the  boundary  terms  and  the  gradient,  that  is 

VL  =  VI  +  [  ^%r7WE  +  [  (42) 

7s  oxj  7S  dxi 

In  the  particular-  case  under  examination  the  functional  to  be  minimized  is 


£ 


U>1  D  +  U>2 


(L  —  L*)2 
2 


(43) 


where  D  is  the  drag  (to  be  minimized),  L  is  the  lift,  L*  is  the  desired  lift  and  the  cy  are  weights.  Consequently, 
the  gradient  expression  is: 


SC  =  cu1J^-(K-%^5y  +  ?^6y)t?dZ+ 

(L  -  L*)J^  -  (~%^Sy  +  ^ Sy )  dS+ 

JtiT%SydZ  +  Jr^SydE 
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5  Use  of  approximate  fitness  evaluators 

Optimization  techniques  based  on  evolutionary  computing  arc  very  attractive  in  terms  of  robustness  and  global 
extremum  location  capabilities,  but  these  favorable  characteristics  arc  obtained  to  the  expense  of  a  rather  large 
number  of  fitness  evaluations.  This  is  a  problem  in  many  engineering  applications,  where  fitness  evaluation 
often  requires  a  substantial  amount  of  computational  resources.  This  makes  attractive  the  use  of  approximate 
fitness  evaluators,  with  lower  computational  requirements,  whenever  it  is  possible.  In  aerodynamic  optimiza¬ 
tion  design,  for  example,  various  approaches  have  been  experimented,  such  as  the  use  of  Neural  Network-based 
interpolators  to  reduce  the  number  of  true  flow  field  evaluations  [12].  The  mixed  use  of  solvers  with  different 
levels  of  approximation  in  a  hierarchical  organization  has  also  been  object  of  investigation.  Both  procedures 
have  advantages  and  limits.  Neural  Networks,  for  example,  arc  able  to  improve  the  quality  of  their  approxima¬ 
tion  when  the  number  of  exact  evaluations  available  increase.  Their  performance,  therefore,  may  improve  in 
real  time  during  the  optimization  process.  On  the  other  hand,  their  more  serious  drawback  is  the  fast  decrease  in 
approximation  fidelity  when  the  parameter  space  dimension  increases  [13,  14,  15,  16].  Low  fidelity  solvers  (e.g. 
Euler+boundary  layer  versus  Navier-Stokes),  do  not  have  this  problem,  but  it  is  not  always  easy  to  understand 
the  limits  of  the  approximate  model.  In  both  cases,  the  problem  that  arises  is  the  difficulty  of  devising  a  good 
and  sound  strategy  to  mix  and  interleave  high-fidelity  and  low-fidelity  fitness  approximations. 

5.1  Problem  analysis 

A  simple  analysis  of  the  implications  of  mixing  fitness  evaluators  with  different  levels  of  precision  within  an 
(evolutionary)  optimization  process  is  here  carried  out.  This  analysis  is  made  using  concepts  of  multi-objective 
optimization  [17,  18,  19],  like  Pareto  fronts  and  dominance  criteria,  even  if  the  class  of  problems  considered 
is  strictly  single-objective.  In  fact,  the  relations  between  the  exact  function  /  and  its  approximated  model  g  is 
analyzed  with  the  following  two-objective  problem: 

{min  /  (x) ,  g(x,X0) 

x  £  X  (44) 

X0  €  p(X) 

where  f(x)  £  +  is  the  function  defined  in  a  domain  X  C  W'  for  which  we  arc  interested  in  finding  the  global 
extrema,  and  g(x,  Xq)  £  5ft  is  defined  in  the  domain  X*  =  X  U  p(X)  ;  p(X)  is  the  power  set  of  X  and 
X()  £  p(X)  represents  a  sampling  of  the  domain  X.  In  other  terms  Xq  can  be  thought  as  the  training  set  of  g, 
e.g.  the  set  of  point  in  which  /  (and  possibly  its  derivatives)  are  known  and  their  values  are  used  to  obtain  g. 
Let  us  have,  for  example,  the  following  function: 

f(x)  =  [ x 2  —  10cos(27rx)  +  10]  a;  £  [—1,  3]  (45) 

andlet  g  (x)  a  polynomial  of  degree  10  that  minimizes  if(xi)  ~  d(xi)}2  to  best  fit  f(x)  in  the  least  squares 
sense  for  a  given  training  set  Xq  (see  Figure  2). 

The  Pareto  Front  related  to  problem  (44),  reported  here  in  Figure  3  for  function  45  and  its  approximator,  is 
useful  to  understand  the  relationship  between  /  and  its  approximating  function  g. 

The  training  set  Xq  is,  for  the  moment,  supposed  fixed  and  unchanged  during  the  optimization  process.  In 
other  words,  Xq  fixed  means  that  no  learning  procedure  is  involved  in  the  approximation  process;  there  are, 
instead,  two  different  models,  from  the  beginning  of  the  process,  the  exact  and  the  approximate  one,  that  can  be 
used  to  search  the  optimum. 

When  using  an  approximated  model  for  fitness  evaluation,  there  arc  two  things  that  have  to  be  taken  into 
account: 

•  Shift  of  extremum  points  between  the  two  models; 

•  Value  shift  between  exact  and  approximate  functions. 
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Figure  2:  f(x )  =  [ x 2  —  10  cos(27tcc)  +  10]  and  g  =  best  fitting  polynomial  of  degree  10. 


Figure  3:  Pareto  front  of  problem  44  with  /  specified  in  eq.  45  and  g  10th  degree  polynomial  approximator. 


The  lesser  arc  these  shifts  in  the  extremum  points,  the  better  and  more  reliable  are  optimization  results  obtained 
using  the  approximate  model.  Furthermore  in  these  conditions  it  is  safe  to  interleave  exact  and  approximate 
evaluations. 

A  generic  Pareto  front  of  a  problem  of  type  (1)  is  reported  in  Figure  4.  Four  different  zones  can  be  found, 
that  are  characterized  by  different  behavior  of  g  and  /  in  moving  towards  their  global  optimum  points.  In  zone 

(I) ,  a  solution  that  dominates  the  initial  one  is  surely  closer  to  both  /  and  g  global  minimum  points.  In  zone 

(II) ,  a  solution  that  dominates  the  initial  one  is  closer  to  the  global  minimum  of  /  and  farther  from  g  global 
minimum.  Conversely,  zone  (IV)  has  an  opposite  behavior,  as  a  dominating  solution  leads  close  to  g  and  far 
from  /  global  minimum  points.  In  zone  (III),  finally,  a  solution  that  dominates  the  initial  one  is  farther  from 
both  global  minimum  points. 
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Given  this  framework,  we  want  check  if  there  is  a  particular  interleaving  strategy  that  minimizes  the  number 
of  exact  function  evaluations.  The  obvious  limit  of  this  approach  is  that  in  a  real  optimization  run  the  Pareto 
front  is  unknown  and  only  a  rough  estimate  of  the  collocation  of  a  particular  solution  in  one  of  the  four  zones  is 
possible. 

At  the  beginning  of  the  optimization  process  the  candidate  solution  (or  population)  is  likely  to  be  in  zone 
(I).  The  approximated  function  g  can  be  safely  used  here,  and  only  occasional  control  is  needed  on  /  values  to 
avoid  the  generation  of  too  many  non  dominated  solutions  (relatively  to  /). 

Zones  (II)  and  (III)  are  difficult  to  detect  if  g  global  minimum  value  is  unknown.  In  principle,  use  of  g 
without  check  on  /  should  be  avoided  in  zone  (II),  as  it  may  lead  to  unsatisfactory  results  in  terms  of  /.  In  zone 
(III)  both  the  use  of  g  alone  and  of  g  and  /  interleaved  is  detrimental,  as  it  would  lead  to  a  point  in  the  Pareto 
front  that  would  be  likely  to  be  far  from  the  global  optimum  of  /. 

Zone  (IV)  is  similar  to  zone  (I)  from  the  point  of  view  of  searching  for  the  global  optimum  of  /.  However, 
here  interleaving  could  be  detrimental  for  reaching  the  global  optimum  of  g. 


Figure  4:  Generic  Pareto  front  related  to  problem  44. 


5.2  Possible  approaches  to  mixing  exact  and  approximated  evaluations 

From  a  practical  standpoint  the  above  considerations  suggest  the  following  search  strategy: 

1 .  Search  the  global  optimum  of  g  first,  eventually  checking  in  very  few  points  if  the  values  of  /  do  not  get 
worse  too  much. 

2.  Search  the  global  optimum  of  /,  eventually  using  g  as  second  objective  if  there  is  interest  in  finding  so¬ 
lutions  belonging  to  the  f-g  front. 

When  Xq  changes  during  the  optimization  process,  the  approximation  g  can  increase  its  precision  when 
more  evaluation  of  /  are  available.  Therefore  the  Pareto  front  tends  to  decrease  in  size,  and  after  some  time  it 
can  reduce  to  a  point. 

Devising  a  good  interleaving  strategy  is  here  more  difficult,  however  the  previous  strategy,  adapted  as  fol¬ 
lows,  should  have  an  acceptably  good  behavior: 
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1.  /is  sampled  in  some  points  and  the  resulting  values  (training  set)  arc  used  to  find  a  first  expression  of  g. 

2.  The  global  optimum  of  g  is  searched. 

3.  A  short  two-objective  ( f-g )  run  is  performed  and  the  computed  values  of  /  arc  added  to  the  training  set. 

4.  A  new  expression  of  g  is  found. 

5.  The  process  is  repeated  starting  from  step  2  until  a  satisfactory  value  of  /  is  found. 

5.3  Analysis  of  an  example  problem 

The  RAE  2822  airfoil  is  assigned  as  starting  configuration,  and  the  optimization  goal  is  the  reduction  of  the  drag 
coefficient  Cd  at  given  mach  (M  =  0.78)  and  lift  coefficient  c/  =  0.60.  The  maximum  airfoil  thickness  has  to 
remain  unchanged. 

A  two-objective  optimization  run  is  carried  out  using  two  different  flow  solvers: 

1.  Euler  +  Boundary  layer,  interactive  (high  fidelity  solver,  objective  /). 

2.  Full  potential,  non-conservative  formulation  (low  fidelity  solver,  objective  g). 

It  is  important  to  observe  that  these  flow  solvers  represent  the  real  flow  at  a  very  different  approximation  level, 
and  that  a  real  optimization  application  would  benefit  more  from  the  use  of  more  closely  related  flow  solver 
models,  such  as  Navier-Stokes  and  Euler  +  Boundary  layer.  The  only  aim  of  this  run  is  the  analysis  of  the  rela¬ 
tionships  between  different  flow  models. 

The  aerofunctions  obtained  from  the  airfoils  reported  in  figure  5  were  used  to  modify  the  airfoil  shape. 


airfoil  #1  airfoil  #2  airfoil  #3  airfoil  #4 


airfoil  #5  airfoil  #6  airfoil  #7  airfoil  #8 


airfoil  #9  airfoil  #10  airfoil  #11  airfoil  #12 

Figure  5:  Airfoils  used  to  build  the  modification  functions. 


The  GA  parameters  are: 

•  bit  resolution:  10 

•  selection  method:  two-step  random  walk  with  elitism 

•  bit-mutation  with  4%  rate 

•  population  size:  20 


generations:  20 
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Figure  6  reports  the  Pareto  front  obtained  in  the  two-objective  run.  As  it  can  be  observed,  the  Pareto  front 
shows  a  sensible  difference  between  the  two  flow  models.  These  solvers  should,  therefore,  not  be  used  in  an 
interleaved  fashion  unless  a  viscous  correction  is  introduced  in  the  full  potential  one.  However,  a  hierarchical 
use  of  these  solver  could  be  useful  to  discard  very  bad  solution  such  as  those  with  a  strong  shock  wave. 

The  above  example  may  help  in  tracing  the  guidelines  for  the  implementation  of  the  hierarchical  approach 
in  an  optimizer  for  viscous  flows. 

The  optimization  could  be  made  in  two  different  steps:  a  pre-optimization  phase  aimed  at  the  reduction  of 
the  search  space  complexity,  and  a  subsequent  optimization  phase  that  explores  this  reduced  space  using  a  high 
fidelity  fitness  evaluator. 

The  number  of  active  variables,  their  variation  ranges  and  a  new  starting  configuration  are  a  possible  output 
of  the  first  phase.  In  a  more  sophisticated  approach,  using  for  example  a  multi-objective  GA,  the  pre-optimiza¬ 
tion  can  be  used  to  generate  an  aerofunction  basis  with  a  reduced  number  of  elements. 

In  the  second  phase,  the  modification  of  the  geometry  could  be  driven  by  a  very  simple  evolution  strategy, 
such  as  a  (1,1)-ES  [20,  21]  that,  due  reduction  of  the  search  space  complexity,  should  give  an  improvement  of 
the  objective  function  in  an  acceptable  number  of  attempts. 


6  Conclusions 

It  has  been  shown  how  computer  algebra  techniques  can  be  used  to  obtain  the  adjoint  equations  and  pertinent 
boundary  conditions  of  the  Navier-Stokes  equations.  The  obtained  equations  can  be  directly  used  for  the  adjoint 
equations  solver  coding. 


full  potential:  cdw=0. 000281 
Euler  +  BL:  cdw=0. 000006 


Figure  6:  Euler  +  Boundary  Layer  —  Full  Potential  comparison. 
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MAPLE  was  used  as  computer  algebra  language,  complemented  with  an  original  code  taking  advantage  of 
public  domain  libraries  for  tensor  algebra  and  function  variation  computation. 

Further  development  will  include  derivation  of  the  adjoint  equations  of  the  Reynolds  averaged  Navier-Stokes 
equations  with  n-e  turbulence  modeling. 

The  relationship  between  an  objective  function  and  an  approximated  form  has  been  investigated  using  a 
Pareto  front  to  show  the  trade-off  between  exact  and  expensive  objective  function  evaluations  and  inexpensive 
but  inexact  approximations.  Approximating  the  objective  function  values,  if  properly  done,  shortens  the  time 
in  time-consuming  design  problems  like  those  typical  of  CFD  applications.  Evolutionary  algorithms  normally 
require  more  objective  function  evaluations  than  other  methods  and,  therefore,  arc  very  likely  to  take  advantage 
by  using  such  approximation  techniques. 
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A  Example  of  adjoint  NS  derivation  using  MAPLE 

For  the  sake  of  clarity  and  conciseness,  the  output  of  MAPLE  has  been  slightly  modified  such  that  indexed  vari¬ 
ables  and  derivatived  such  as  xl  and  v‘oX2  appeal-  as  x\  and  vj,xr  Furthermore,  we  will  work  only  on  a  piece  of 
the  functional  20,  namely  the  scalar  products  between  the  adjoint  variables  and  the  state  equations.  Therefore, 
in  the  output  will  not  appeal-  the  terms  coming  from  I.  Hence  we  have  the  following  expanded  form: 


LG1  :  = 


V  ( Px i  Vi  +  p  vlxi  +  px&  v2  +  p  vZx2  +  pXs  v3+  p  v3xs)  +  A i  (pXl  Vi' 


4  1  1 

+  2 pvi  viXl  +pXl  -  -nviXuXl  -  -pv2xi,X2  -  ~pv3xuxs  +  pX2  v2  vi  +pv2xz  vi 
+  pv2  Vlx2  -  pviX2jX2  +  pXs  v3  Vi  +  pv3xs  Vi  +  pv3vlxs  -  pvlxSiX3)  +  X2 

,  1  9 

( pxi  V2  Vi  +  pv2xi  Vi  +  pv2viXl  -  pv2xuXl  -  -pviXuX2  +  pX2  v2z  +  2  pv2v2xe 

4  1 

+  Pxz  —  2  P  V%x2,  x2  —  2  P  v3  x% ,  x$  +  Pxs  V3  V2  +  P  v3xs  V2  +  p  V3  V2xs  —  p  V2xg^  X3)  + 
^s{Px!  v3  Vi  +  p  v3xi  Vi  +  p  V3  ViXl  —  p  v3xi  j  Xl  —  —  p  Vi  Xl )  Xs  +  pX2  v3  v2 

1  9 

+  P  v3x2  V2  +  p  v3  V2xz  -  p  v3x2jX2  -  -  p  v2xSt  xs  +  pxs  v3z  +  2  p  v3  V3xs  +  pxs 

4  14 

-  3  PV3xgyX3)  +  C(-g  v3pvixuxs  -  -  V3xs2  p-  V2xs2  p  +  pX3  v3  H  +  p  V3xs  H 

T  P  v3  HXs  —  k  TX3j  xs  —  Vi  X3  p  —  k  TXl )  Xl  —  —  v2  p  Vi  Xl )  X2  —  k  TX2t  Xg  —  —  v2  p  v3xgt  Xg 

4  1  4 

Vl  /i  v  1x2,  X2  2  M  ,X2  ^3  M  ^3  Xl  ,Xi  2  Vl  ^  ix2  2  V3  ft  ^ 3  Xs  ,  Xs 

1  4 

-  2  v3pv2x2jXs  -  Vl  pviX3tX3  -  v3pv3X2}X2  -  -  ViXl2  p  +  pvi  HX1  +pviXl  H 

4 

+  Pxi  vi  H  -  v3x2  p  -  v2pv2x3tXs  -2v2xi  pviX2  -  v2xi2p  +  -  viXl  pv3xg 

4  4  1 

“I-  2  xi  M  ^2  X2  ^2  M  ^2  xi ,  xi  2  xi  M  xs  2  M  Vl  xi ,  xi  7^^!  [l  V3  Xl  ?  X3 

4  4 

-  3  v2x2  p  -  vlx2  p  +  pv2HXz+  p  V2x2  H  +  pX2v2H  -  v3x2  p  +  -  v2xg  p  V3xs 


-2v3x2  pv2xs)dxidx2dx3  +  J J  £1  vi  +  i2  v2  +  £3  v3  +  r@(T)dsi  ds2 
Now  the  perfect  gas  model  is  introduced  in  order  to  close  the  equation  system: 

>  pexp  :=  p  =  -l/2*rho* (gamma-1) * (-2 *H+vl ~ 2+v2 ~ 2+v3 ~ 2 ) /gamma; 


pexp  :=  p  =  — - 


1  p  (7  -  1)  (-2  H  +  vi 2  +  v22  +  v3 2) 


7 


>  Texp  :=  T  =  -1/2* (gamma-1) * (-2 *H+vl ~ 2+v2 ~ 2+v3 ~ 2 ) / (gamma*R) ; 

1  (7  -  1)  (-2H  +  vi2  +  v22  +  v32) 


Texp  :=  T  = 


7  R 


>  LG1 : =subs ({pexp, Texp},  LG1 ) : 

The  prescribed  boundary  temperature  condition  is  imposed: 

>  Theta  :=  (w)  ->  w -  Tw(xl,x2,x3,sl,s2); 

0  :=  w  — >•  w  —  Tw (xl,  x2,  x3,  si,  s2) 

The  variation  with  respect  to  p  gives  the  first  adjoint  equation: 

>  ad j_drho : =compute_ad j_eq (LG1 ,  rho,  drho,  2,  0,  true); 
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rjf  1  V3  A 3xg  \  \  H  ^3xg  .  H  ^2xg  \lXl  H 

~v2  H  QX2  -  - - -  -  v2  V3  X3x2  -  v2  V3  X2xs  H - -  H - -  -  Vxs  v3  H - 

27  7  7  7 

-  Aia;j  H  -r)Xl  v !  -  H  X3xs  -  pX2  v2  +  l;  vi2\3xs  -  ^  vi2 \iXl  +  ^  vi2  X2xg 

+  2  ’u52  ^2xg  +  ^  A^  xi  ^3“  +  2  +  2  ^5:rs  ~  2  1,22  ^2x2  _  2  1,52  ^5xs  _  H 

rj,  1  Aljj  Wj2  ,  1  V!2X3x3  1  V32  X2x2 

~V!H  QX1  -  - V 1  V3  Xlx  -  - -  - -  -  V 1  V3  X3xi 

27  2  7  2  7 

1  v/2  X2x  1  V!2X1x  1  v22X3x 

-  77  - -  -  77 - (xsV3H  -  -  - -  -  V!  V2  X !  X2  ~  Vi  V2  X2xi 

272727 


.2  1  „,_2 

2  7 


1  A iXl  v32  1  n22  X2xs 


2  7 

Now  the  adjoint  equations  coming  from  the  variations  with  respect  to  v\,  v2,  '('3,: 


>  ad j_dvl : =compute_ad j_eq (LG1 ,  vl,  dvl,  2,  0,  true); 


5  4 

~P  H  Cxi  ~h  g  P  V2  xi  Cxg  ~  P  Cxg,  xg  vl  —  g  P  V1  Cxi  ,xi  —  p  ^lx2,xg  ~  P  v3  ^lxg 


+  pvi  X2x2  —  pv3  X3xi  — 


P  vi  A 1  Xl  p  vi  X2x2  p  vi  X3xg  k  vi  C Xl .  Xl 


7 


7 


7 


+ 


R 


kvi  C Xg,Xg  ,  kCxs,X$  V1  5  5 

H - ^ ^ - Vxip-^p  V3xs  (X1  -  ~  p  Cxi  v2xe 

—  2  P  v3  Cxi  ,xs  4“  2  P  v3 xi  Cxg  g  P  A;  xi ,  xi  4“  P  vi  X3 Xg  —  p  X2 Xl  )  Xg  P  V2  A;  X2 


1 


1 


o  P  V2  Cxi  ,  Xg  nP^3xi,Xg  PV1  XlXl  p  Vl  yHX2  X2  pXlXgXs 


k  Cxg,  Xg  Vl 

7  R 


k  Vl  Cxi  ,X1  k  Vl  Cxg,  Xg 


7  R 


7  R 


—  pv2X2 


Xl 


>  ad j_dv2 : =compute_ad j_eq (LG1 ,  v2,  dv2,  2,  0,  true); 

1  5  5 

—  p  Cxi  ,X1  V2  —  -  p  V3  Cxg,  xg  +  p  V2  X3xs  +  P  v2  XiXl  +  -  p  V3xg  Cxg  +  ^  P  Vl  X£  Cxi 

,  pv2Xix  pv2X3x  pv2X2x  4 

-PX2xi,Xi - - - - - - Pxgp-  ^pX2xg,Xg 

k  C  Xl ,  Xl  V2  k  V2  Cxg  ,Xg  k  V2  Cxg  ,xg^-,  4 

quR  7 R  qi?  3M  3x2,xs  ~  ^Pv2  Cx2,x2 

5  1 

—  p  V3  X2xs  —  —  p,  Cxg  Vl  Xl  —  p  V2  X2x2  —  pv2  Cxg,  Xg  ~  P  H  Cxg  ~  ^  p  V1  Cxi  ,  Xg 

\  \  1  \  \  .  kv2Cxg,Xg  .  k  V2  Cxg ,  Xg 

—  pv3  X3x2  —  pvi  X2xi  —  -  pXiXuX2  —  pvi  XiX2  H  — - 1  - 

k  C.xi  Xi  ^2  \  5  . 

4  P  ^2xg,Xg  ~  g  P  V3 xg  Cxg 

>  ad j_dv3 : =compute_ad j_eq (LG1 ,  v3,  dv3,  2,  0,  true); 
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2  P  Cxi  ,xs  Vl  g  P  "^2 X2,xs  P  V 2  X2 xs  P  ^3  X 3 xs  P  X}  x3  P  H  Cxg  g  P  ^  1  xi ,  xg 

,  ,  PV3^2xz  Pv3^1Xl  PV3X3 

P^3xi,Xl  P^3xg,X2  Vxg  P  ^ 

k  Cx2  ,X2  V3  k  C Xl  ,  Xi  T?  k  V3  Cxg  ,Xg  .  k  C,X1  ,  xi  L? 


xg 


+ 


3P^3xg,Xg  iR 

k  Cx2  ,X2  ^3  k  V3  Cxg ,  Xg 


7  R 


+ 


R 


+ 


R 


R 


5  5 

P  ^2  X3  X2  “1“  g  P  Xg  Cxi  g  P  ^x3  ^  X1 


5  5 

—  P  Vi  \sXl  +  p  V3  \iXl  —  -  pCxg  V2X2  +  2  P  ^Xs  V2xs  P  V3  ^2xe  ~  V  I  xs  v3 

4  1 

—  2  PV3  C Xg  ,Xg~P  Cxi  ,X1  v3  ~  2  P  V2  Cx2 ,  Xg 

Finally  the  variation  with  respect  to  H  gives: 

>  ad j_dH : =compute_ad j_eq (LG1 ,  H,  dH,  2,  0,  true); 


P  X3  xg  + 
0 

+ 


k  C 


Xl  ,  Xl 


7  R 

P  1  Xl  k  CtXl  ,  Xj 


P  T 2  Cx2  P  ^2x2  P  Vl  C Xl  P  L?  Cxg  4" 

k  C Xg ,  Xg  .  P  X 2x2  ^  Cxg  ,  X2 


k  Cxg, 


xg 


k  C 


X2,X2 


P  X 1  xi  + 


7  R 

P  ^  3  xg 


7  R 


7  R  R  7  R  *  7 

These  equations  can  then  be  combined  in  order  to  obtain  the  variation  related  to  the  conservative  variables 

P ,  pvi ,  pE. 

The  same  process  is  followed  to  obtain  the  boundary  conditions.  However,  in  this  case  there  is  a  further 
complication,  because  the  equations  have  terms  containing  the  variations  (dp,  Svi,  5H).  See  sub-section  4.3, 
and  eq.  38  for  a  proper  treatment  of  such  terms. 


B  Formulas  for  functional  variations 

B.l  Variation  of  surface  integrals 

In  the  following  section  the  variation  of  a  surface  integral  (see  eq.  42)  with  respect  to  surface  changes  is  explicitly 
derived,  in  a  form  useful  for  computer  algebra.  Let 


S  =  x(u,v)i  + y(u,v)]  +  z(u,v)k 

be  a  parametric  representation  of  the  integration  surface.  The  surface  integral 

f(x,  y,  z )  da 


(46) 


can  be  expressed  as  a  double  integral  as  follows: 

pb  pS( v) 


a 


f(x(u,  v ),  y(u,  v),z(u,  v))W (u,  v)  du  dv, 


where 


with 


Let 


W(u,v)  =  yfj[+j[+Ji 


Jl  = 


d{y,  z) 


d{u,  v ) 


■h  = 


d(z,  x) 


d{u,  v ) 


■h  = 


d(x,y) 


d{u,  v ) 


S'  =  [: x(u,v )  +  hi(u,v,a)]  i+  [y(u,v)  +  h-2(u,  v,  a)]  j  +  [z(u,v)  +  h3(u,v,  a)]  k 
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an  arbitrary  variation  of  S  controlled  by  a  generic  parameter  vector  a.  Now  we  can  write: 


rb  rS(v) 

F(  a)  =  /  /  f(x  +  hi,y  +  h2,z  +  h3)W(u,v,a)dudv. 

J  a  J  j(v) 


The  integration  boundaries  arc  unchanged  and  independent  on  a.  Therefore  only  the  variation  of  the  integrand 
has  to  be  considered  when  computing  the  derivative  with  respect  to  a  generic  element  a  of  a: 


dF{  a) 
da 


fb  fs 0)  <9h(a) 

/  /  V/- — - — -W(u,  v,  a)  dudv+ 

'a  J-y(v)  Oa 

r  r  *>*>■ 


a  Jry( v ) 


If  /  depends  explicitly  on  a,  then  a  further  term  has  to  be  considered,  and  the  above  equation  becomes: 


dF(  a) 
da 


fb  f5F)  <9h(a) 

/  /  V/- — - — -W(u,  v,  a)  dudv+ 

la  va 

[b  [s^  dW(u,v,  a)  ,  ,  fb  [5{v)dfTXr,  ,  J  J 

/  /  / - - - -dudv+  /  /  —W(u,v,a)dudv. 


a  J  7(11) 


a  J  7(f) 


Considering  that 


dW(u,  v,a)  Ji  Jia  +  J2 J2a  +  J3J3C,  JiJia  +  J'2J2a  +  J3J30 


'J\  +  J22  +  4 


Jf  +  Jf  +  Jf 


W  =  HW 


we  finally  have: 


dF{  a) 
da 


as{-v)  dh.(a) 

V/  •  — - — -FF(u,  n,  a)  dudv+ 

» 

rb  l*S(v)  rb  rS(v)  q  j 

/  /  fF[(u,v,a)W(u,v,a)dudv  +  /  /  —  W(u,v,a)dudv. 


a  J  7(1;) 


a  J  7(1;) 


It  can  be  shown  (see  next  sub-section)  that  H  does  not  depend  on  the  particular-  parametric  representation  chosen, 
and  it  is  therefore  possible  to  write,  finally: 


dF(  a) 
da 


V/  •  —  dcr  +  /  /  fH  da  +  /  /  fada. 


In  terms  of  variation  of  h  the  above  formula  can  be  expressed  as: 


SF=  Vf -6h  da  +  //  fSnhda+  8fhda 


J  j  S  J  j  ji . 

qj  _ _  L'h 

3~  Jf  +  Jl  +  Jl 

Finally,  in  terms  of  variation  of  a  generic  function  g,  we  have: 


8F  =  jj  V/  •  Sg  da  +  JJ  fU  ■  5Ug  da  +  JJ  ^  Sg  da 

s  s  s 
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B.2  Invariance  of  H  with  the  parametric  surface  representation 

To  demonstrate  the  invariance  of  TL  j  it  is  sufficient  to  show  that 


Oj(s,  t) 

Ji(u,v ) 


=  usvt  -  utvs 


and 


dJj(s,t,  h) 
dhj 

dJj(u,v,  h) 
dhi 


=  usvt  -  utvs 


This  can  be  achieved  with  simple  algebraic  passages: 

9 


d{y,  Z ) 

d(u,  v ) 


d^y(uiv)  aT,y(u’v) 

z(u,v)  Sr,z(u’v) 


d_ 

L  9m 


d(u,  v ) 
d(s,t) 


Ts  u(s>  t) 


m  u(s>  t) 


£v(s,f)  4v(*,t) 


d{y,  z)  =  d{y,  z)  _  Q(n,n) 
0(s,  f)  0(n,  v)  0(s,  t ) 


liyfu  +  l;y^v  l^ylu  +  l;y|v 


d  „  d 


d  „  9 


9m  z  9s  u  +  9m  z  9s  v 


9^9 


a  „  a 


du  z  at u  +  9m  z  at v 


Ji(s,t )  0  0 

Ji(u,v )  Os  Of 


0  0 
Os  U  Of  V 


OJi(n,  n) 
dhi 


02 


On  dhj  dv 


d  d 
y^r  z  + 


02 


02 


On  On  dhi 


■  z  — 


On  dhj  du 


d  d 
Yttz-  ^y 


02 


On  On  dhi 


d{dhj’z ) 

0(n,  n) 


92 


9m 


92 


dudhi  y  dvdhi  Y 


ikz 


o  (y.  fj 

0(n,  n) 


&y 

92  . 
dudhj 


*y 

92  , 
dvdhj 


dJi(s,  f,  h) 


dhi 


aiyhpz)  9(m,m) 
9(m,m)  d(s,t) 


+ 


d(y,Zhj)  9(m,m) 
9(m,m)  d(s,t) 


92 


9  „  9 


a2  VAUAZAV+  a2  vlvlzlu _ 9i_  viulzlv _ v  JL  v  JL  z 

dudhj  y  9s  u  9m  z  dt  V  ^  9m 9/i,  Y  9s  v  9m  z  dt  Ll  dudhj  Y  dt  Ll  9m  z  9s  v  9m 9/i,  Y  at  V  9m  z 

9,,9..  92  „  9  „  ,  9  „  9  „  92  „  9  „  9  „  9  „  92  „  9  „  9  „  9  „  92 


911  y  95  u  99M-  z  m  v  +  &  y  95  v 


dhj  dt 


z  *  u  -  y  ft  u  917M7  z  93  v  -  KJ  y  m  v  dudhi 


dJj(s,t,  h) 

dhj  d  d  d  d 

dJi(u,  v,  h)  Os  V  Of  U  Os  U  Of  V 

dhj 
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B.3  Volume  integral  variation 


The  volume  integral  variation  is  obtained  making  explicit  reference  to  the  reduction  formulas  of  a  volume  in¬ 
tegral  into  a  triple  integral,  and  applying  in  a  straightforward  way  the  formula  for  the  derivative  of  a  definite 
integral: 

rv(t)  Mt)  g  d  d 

f(x,t)dx=  I  -f(x,t)dx  +  f(v(t),t)—v(t)-f(u(t),t)—u(t) 


lu(t) 


lu(t) 


dt 


dt 


which,  applied  to  a  triple  integral,  gives: 

Q  rb(t)  rS{x,t)  M(x,y,t) 
dt  J a(t)  J<p(x,y,t) 


f(x,  y,  z,  t)dzdydx  = 


(■b(t)  rS(x,t)  rij)(x,y,t)  Q 


'  a(t)  7 (x,t)  J(j)(x,y,t) 


g^f{x,y,z,t)dzdy  dx 


+ 


V,  t )  f(x,  y,  tp(x,  y,  t),t)dydx 


cb(t)  rS(x,t)  /  Q 
i(t )  \dt 

fb(t)  r8{x,t)  /  Q  \ 

/  t)  /(*>  Vi  Vi  t),t)dyda 

i(t)  JMx,t)  \at  / 


+ 


' a(t) 

C(i‘M 


ri/j(x,8(x,t),t) 
/  (j)(x,5(x,t),t) 


f(x,  6(x,  t),  z,  t)dzdx 


I  a(t ) 


n'i/j(x,,y(x,t),t) 


I  (j)(x^(x,t),t) 


f(x,  i(x,  t ),  z,  t)dzdx 


+  \jtb{t) 


rS(b(t),t)  M(b(t),y,t) 

'7 (b(t),t)  J<t>(b(t),y,t) 


f(b(t),y,z,t)dzdy 


a(t ) 


nS(a(t),t)  1-4 >(a(t),y,t) 
'7 (a(t),t)  J <p(a(t) ,y ,t) 


f(a(t),y,z,t)dzdy 


Introducing  the  deformation  field  h(x,  y,  z),  the  boundary  function  ?/;  and  (j)  can  be  written  as: 

^(x,  y,  t )  =  ip0(x,  y)  +  h3(x,  y,  i/>0(x,  y),  t ) 

4>(x,  y,  t )  =  4>q(x,  y)  +  h3(x,  y,  (j)0{x,  y),  t ) 

This  allows  the  identity  below: 

Mt)  rS(x,t)  /  q  \ 

/  /  [wL‘llJ(x,y,t))f(x,y,ilj(x,y,t),t)dydx  + 

J a(t)  J'y{x,t)  \Ot  J 

Mt)  rS(x,t)  /  Q  \ 

~  /  -37 8(x,y,t))  f(x,y,(j)(x,y,t),t)dydx  = 

J a(t)  Jl{x,t)  \Ot  J 

Mt)  rS(x,t) 

/  /  Di(h3)(x,y,ip0(x,y),t)f(x,y,ip(x,y,t),t)dydx  + 

J a(t)  Jj(x,t) 

rb(t)  rS(x,t)  f  f  QU 


' a(t)  Jry(x,t) 


Di(h3)(x,y,(j)o(x,y),t)f(x,y,(j>(x,y,t),t)dydx  =  l/—fdydx 


+FT 


Reasoning  in  a  similar  way  for  the  remaining  couples  of  integrals,  we  can  finally  write: 


d_ 

dt 


fdV  = 


d£ 

dt 


dV  + 


dhi 

~dt 


fdy  dz 


dh2 

~dt 


fdz  dx  - 


dh:> 

dt 


fdx  dy 


+FT 


+FT 


+FT 
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That,  in  terms  of  surface  integrals,  becomes: 


d_ 

dt 


fdV  = 

T  T 


FT 


Finally,  in  terms  of  variation  of  a  generic  function  g ,  we  have 

5F  = 


j-gdgdV  +  II  f^.ndgda 
T  FT 


B.4  Application  of  divergence  theorem  to  terms  involving  crossed  derivatives 

The  divergence  theorem 

III  V  •  b  (in  =  II  b  •  n  dT, 

n  s 

can  be  applied  in  two  different  ways  to  transform  the  following  volume  integral 


(48) 


dfl 


with  v  a  generic  scalar  function  continuous  in  0  with  its  first-  and  second-order  partial  derivatives,  such  that 
vxy  =  VyX .  It  is,  indeed,  possible  to  write  vxy  either  as  V  •  (0,  vx,  0)  or  as  V  •  (vy,  0,  0).  This  leads  immediately 
to  the  following  identity: 


Txy  r / -1 


vxri2  dTi 


vyrii  dTi 


n 


s 


s 


with  n  =  (m,  n2,  n3). 

This  identity  can  be  used  to  transform  some  adjoint  boundary  conditions.  This  operation,  indeed,  can  be 
useful  to  restore  symmetry  in  the  adjoint  boundary  conditions  automatically  derived. 
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